Functional MRI-derived priors for solving the EEG/MEG inverse problem 
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1 Introduction 

Because of their excellent temporal accuracy (of the 
order of 1 ms), electroencephalography (EEG) and 
magnetoencephalography (MEG) provide the most 
relevant data for studying the temporal dynamics of 
brain activity. However, difficulties arise when trying 
to localize the electromagnetic sources of this activity 
from EEG/MEG scalp recordings. This mathematical 
inverse problem is indeed ill-posed and largely un¬ 
derdetermined. An efficient way of constraining the 
problem and thereby reducing the solution space is to 
perform a regularization. By taking some anatomical 
and/or functional a priori knowledge into account, the 
regularization process may yield a more consistent lo¬ 
calization of the electromagnetic sources. Anatomi¬ 
cal priors have already been used (e.g., [1]) but only 
few regularization methods (e.g., [2]) have introduced 
functional information so far. 

In this study, we propose a new multimodal ap¬ 
proach for solving the EEG/MEG inverse problem. 
This method involves a distributed source model and 
accounts for anatomo-functional constraints derived 
from functional magnetic resonance imaging (fMRI) 
data. In the following, we briefly describe the source 
model, the regularization procedure and the way func¬ 
tional priors are introduced. In order to assess the 
value of the proposed approach, we then present re¬ 
sults obtained using simulated data. 

2 Method 

The brain activity localization method used in the fol¬ 
lowing is based upon the ST-MAP (Spatio-Temporal 
Maximum A Posteriori) algorithm developed by Bail- 
let and Garnero to solve the inverse problem in 
EEG/MEG [1]. In section 2.1, we introduce the main 
features of this algorithm, namely the distributed 
source model and the regularization process involv¬ 
ing spatial and temporal priors. In section 2.2, we em¬ 
phasize the new developments that have been made in 
ST-MAP in order to account for functional priors. 


2.1 Distributed source model and regular¬ 
ization procedure 

We wish to localize brain activity electromagnetic 
sources from electric and/or magnetic scalp mea¬ 
surements covering a time range of a few hundred 
of milliseconds. In the so-called distributed source 
model, these sources are conventionally modeled as 
current dipoles with given positions (it is assumed 
that dipoles may be located only on the cortex strip 
[3, 4]) and orientations (since it has been shown that 
neuron columns are organized perpendicularly to the 
cortical surface [4]). 

Using segmented high-resolution anatomical MR im¬ 
ages, we initialize the localization process by spread¬ 
ing N ~ 7,000 dipoles (whose position and ori¬ 
entation are fixed) over the interface between white 
and gray matter (Fig. 1). The amplitude of each 
dipole is the only parameter that still remains to be 
determined. The inverse problem then consists in 
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Figure 1: a) Segmented white matter image (axial 
slice), b) Position of dipoles over the cortical surface. 

finding the N amplitude values that best explain the 
EEG/MEG measurements. This can be expressed us¬ 
ing the following linear system: 

M n = GJ n + b n , (1) 

where M n is the vector containing the EEG/MEG 
measured data at time step n, G is the forward prob¬ 
lem operator, J n is the vector containing the nth time 




sample of dipole amplitudes and is an additive 
measurement noise. 

Since the problem is ill-conditioned, a regularization 
procedure is mandatory. However, applying the regu¬ 
larized localization procedure on all the N sources is 
untractable: the number of initial dipoles is too large 
to yield a coherent solution and the computation time 
is unacceptable. We therefore select only N ~ 2,000 
dipoles for solving the inverse problem. This choice 
is guided by analyzing the contribution to the data 
(lead field) of each source considered separately. 

The regularization procedure used for recovering 
electromagnetic sources is then based on a Bayesian 
approach which consists in minimizing the following 
criterion (Eq. 2): 


||M n — GJ n || 2 
^ ^ ^ ' &K P n/i ^n/j) 
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where 3 n / i is the amplitude of the ith dipole at time 
step n, Ci is a spatial neighborhood of the ith dipole, 
A is a regularization parameter (hyperparameter), <fi K 
is a potential function and P^-i is the projector onto 
the space orthogonal to J n _i. 

The main part of the criterion, ||M^ — GJ n || 2 , rep¬ 
resents the data likelihood term. The remaining terms 
are prior terms whose influence depends on the value 
of the hyperparameter A. Two types of constraints can 
be distinguished: 

• ||P^_ 1 J n || 2 represents the temporal prior. This 
term enforces the continuity between the ampli¬ 
tudes of a given dipole estimated at two consec¬ 
utive time steps. 

• E ( E ] represents the spa- 
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tial constraints. The potential function <fi K is a 
Lorentzian function which is applied to the am¬ 
plitude gradient of two dipoles i and j belong¬ 
ing to a common spatial neighborhood Q. The 
behavior of <fi K and the influence of the param¬ 
eter K are illustrated on Figure 2. Minimiz¬ 
ing 4>x(u) corresponds to enforcing the correla¬ 
tion between the amplitudes of two neighboring 
dipoles (u denotes the amplitude gradient). For 
a given pair of dipoles, K is set to a fixed value 
based on anatomical information (relative posi¬ 
tions and orientations of the two dipoles on the 
cortical surface). 



Figure 2: The Lorentzian potential function. 


Regularization is achieved by using the ARTUR itera¬ 
tive algorithm [5], which is initialized using a simple 
minimum-norm solution. Even though about 2,000 
dipoles are initially selected over the cortical surface, 
the whole ST-MAP method may yield a final solution 
with only a few activated dipoles since the least sig¬ 
nificant sources are discarded at each iteration step. 

2.2 Functional priors 

Its ever-improving spatial resolution makes fMRI 
the best candidate for adding functional priors in 
EEG/MEG methods of source localization. It is in¬ 
deed natural to exploit the complementarity of differ¬ 
ent brain functional exploration techniques, namely 
to combine the good spatial resolution of fMRI with 
the excellent temporal resolution of EEG and MEG. 
However, since these techniques have their own phys¬ 
ical characteristics, using fMRI and EEG/MEG in¬ 
formation together remains problematic. Since they 
measure different physiological aspects of brain func¬ 
tion, the location and extent of fMRI activated areas 
might differ from those of the EEG/MEG sources. 
Furthermore, since simultaneous recordings are still 
at their very beginning, EEG and fMRI data are 
usually acquired in two separate sessions. Care 
must therefore be taken to use the same experimen¬ 
tal paradigm and errors may arise when registering 
the data obtained with the different techniques to 
the same coordinate system (e.g., that of the high- 
resolution anatomical MR image of the subject). 

For all these reasons, we chose not to introduce too 
strong fMRI-derived constraints, but rather to slightly 
influence the regularization procedure. As in [1], we 
thus avoided to constrain directly the dipole ampli- 









tudes by modifying the main term of the criterion (Eq. 
2). Instead, we modified the spatial regularization 
term as follows: 

EE ( 1 + 7 WiWj) <p K (J n /i - J n/j) ),( 3) 
i =i ^jeCi ' 

where w h and Wj are probability coefficients assigned 
to dipoles i and j respectively, depending on the local 
activity detected using fMRI, and 7 is a hyperparam¬ 
eter. For a given dipole 7 w t is determined as follows: 

• A cylinder with a fixed volume is defined around 
each dipole. A voxel of the MR image is consid¬ 
ered to belong to the cylinder if its center belongs 
to the cylinder. 

• The weight w h assigned to the ith dipole is then 
defined as follows: 


^ - T/ ’ 

not 

where V^t is the number of activated voxels (de¬ 
tected using fMRI) which belong to the cylin¬ 
der and Vt 0 t is the total number of voxels in the 
cylinder, w, thus ranges from 0 to 1 . 

If both neighboring dipoles i and j lie outside an ac¬ 
tivated region detected with fMRI, wi = wj = 0. If 
the ith dipole lies inside an activated region while the 
jth dipole lies outside an activated region, wi ^ 0 
but Wj = 0. Consequently, in both previous cases 
the spatial constraints term is not modified compared 
to that of Eq. 2, which means that only anatomi¬ 
cal prior knowledge is accounted for. On the other 
hand, if both dipoles i and j lie inside an activated re¬ 
gion, WiWj > 0 and the spatial prior associated to 
these dipoles is reinforced, i.e., the correlation be¬ 
tween their amplitudes is enhanced. This allows a 
weighting of the spatial prior term based on functional 
knowledge. 

3 Application 

3.1 Simulations 

A 3D high-resolution (voxel size: 0.9375 mm x 
0.9375 mm x 1.5 mm) MR image from a healthy vol¬ 
unteer was segmented. The boundary between white 
and gray matter was approximated with small trian¬ 
gles whose vertices gave initial dipoles position for 
the distributed source model [ 6 ]. 


An fMRI activated area was simulated at an arbitrary 
location by delineating a parallelepipedic binary mask 
(inside: activated; outside: unactivated) in the right 
hemisphere, using voxels of the 3D MR image (Fig. 
3a). MEG data were then simulated by locating artifi¬ 
cial current sources on the cortical surface: 4 dipoles 
were located within the fMRI activated area and 2 
dipoles were located on the left hemisphere surface 
(Fig. 3b). The 6 dipoles had a normalized ampli¬ 
tude modulated in time by a sine function. Their con¬ 
tribution (lead field) was calculated by applying the 
forward problem operator. This yielded MEG mea¬ 
surements in a ten-time-sample-wide window whose 
central part (time samples 4 to 6 ) was used as the in¬ 
put data for solving the inverse problem (Fig. 4). 
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Figure 3: a) Simulated fMRI activated area (black 
parallelepiped), b) Simulated dipoles (black circles). 



Figure 4: Simulated MEG data. 


To assess the usefulness of accounting for functional 
priors, results obtained using ST-MAP without ap¬ 
plying functional constraints (method 1 ) were com- 

















pared to those obtained by including functional con¬ 
straints derived from the simulated fMRI activated 
area (method 2 ). 

3.2 Results 

Figure 5 shows the 2 dipoles estimated using method 
1 (A = 10 -5 ) and the single dipole estimated using 
method 2 (A = 1CT 5 , 7 = 10). Method 1 failed to 
localize either set of dipoles. For the same configu¬ 
ration (i.e., same A value), method 2 did not localize 
any dipole in the left hemisphere either but correctly 
localized a dipole within the fMRI activated region, 
near the 4 sources simulated in the right hemisphere. 
These results clearly show that functional constraints 
reinforced the spatial prior term and demonstrate in 
this particular case that functional priors may be the 
decisive information for recovering activated sources. 
Yet, the extent of the right hemisphere activated re- 



Figure 5: Simulated dipoles (black circles) and 
dipoles estimated using method 1 (white squares) and 
method 2 (white circles). 

gion was not fully restored using method 2 (a single 
dipole was localized while 4 were simulated) and the 
left hemisphere activated area was not detected (as 
in method 1). However, the results show that func¬ 
tional constraints derived from another exploration 
technique such as fMRI may yield part of the solu¬ 
tion, whereas methods which do not use such prior 
knowledge are misleading. 

4 Conclusion 

We showed how fMRI-derived priors could be added 
when solving the EEG/MEG inverse problem and 
presented promising preliminary results. However, 
several parameters and processes are involved and 


some issues therefore remain to be addressed: 1 ) the 
method is particularly sensitive to the hyperparame¬ 
ters value, 2 ) modifying the potential function (e.g., 
choosing a convex function) and the preliminary pro¬ 
cess of source selection should improve the benefit of 
using functional priors, 3) selecting only the most sig¬ 
nificant sources might hinder a better recovery of the 
extent of activated areas, 4) other ways of expressing 
functional priors could be considered. Other simu¬ 
lation configurations thus remain to be investigated 
before applying this method to real data. 
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